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Abstract 

Miniature cryogenic coolers operating on the reversed Stirling cycle are used for cooling the detectors of thermal imaging systems. 
The type is a candidate for duty in satellites—a singularly demanding environment, since whatever the cooling power achieved, 
there is pressure for it to be greater; whatever the electrical power consumption, there is demand for it to be less. Stirling cycle 
machines, both prime movers and reversed-cycle, are renowned for the high cost of trial-and-error development. The situation calls 
for improved understanding of the cycle of thermodynamic processes. This paper makes a start by breaking with the tradition of 
arbitrary linearisation of the defining equations, and by integrating the characteristics of the electromagnetic drive circuit with 
simultaneous computation of the dynamic, thermodynamic and fluid flow processes. A specific cooler is defined and simulated 
performance points adduced in evidence of correct functioning of the computations. © 1999 Elsevier Science Ltd. All rights reserved. 
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Notation 

a(l,l), b( 1) coefficients of array - 
etc. 

< 2 , b, aa, bb constants defined in text () 

A ff free-flow cross-sectional area m 2 

A w wetted area m 2 

A x cross-sectional area exposed to gas 

pressure m 2 

C f friction factor-equivalent wall shear stress 

divided by 1/2 pu 2 - 

COP coefficient of performance—ratio heat 

lifted/work in - 
D cylinder bore m 

F l/2C f u 2 /r h m/s 2 

E sol force generated by solenoid N 

F' F Ik V 1/3 

' sol ' sol' Xp V sw 

/ cyclic frequency s _1 

h coefficient of convective heat transfer 

W/m 2 K 

H enthalpy J 

k thermal conductivity W/mK 


* Tel: + 44-01223-332653; fax: + 44-01223-332662; e-mail: 
ajo@eng.cam.ac.uk 


Lr 

length of regenerator packing m 

Aef 

reference length ( = V sw 1/3 ) m 

M 

total working fluid mass kg 

m 

variable mass kg 

m' 

mass rate kg/s 

^MA 

speed parameter or characteristic Mach 
number wVI^NrT - 

N pr 

Prandtl number, ptc p /k - 

N RR 

characteristic Reynolds number 

^ma'^sg 

N re 

local, instantaneous Reynolds number 
Apurjpt = 4r h m'/pA f f - 

n sg 

Stirling number p ret /m/x - 

P 

absolute pressure Pa 

q, q' 

heat, heat rate J, J/s 

P, Q 

abbreviations for terms arising in the 
integrating factor method of solution of a 
first-order differential equation (definitions 
in text) (-) 

Q' 

volume flow rate m 3 /s 

r h 

hydraulic radius—free-flow area/wetted 
perimeter m 

R 

gas constant for specified gas J/kgK 

R e 

Annand’s Reynolds number, IfSpD/ /jl - 

S 

stroke m 
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S 

T 

t 

U 

uo 

u 

u 

Vk, V c , 


V 

v sw 

w 

Y 

1 sep 


ampd 

Y 

1 amp p 

Jd 

5 

l± 

P 

cr, a' 


co 

Iv 


specific entropy J/kgK 
absolute temperature K 
time s 

internal energy J 
step function defined in text - 
bulk velocity, m'/pA m/s 
dimensionless u defined u/wL ref - 
V buff volume excursion of 

expansion/compression/buffer space m 3 
swept volume m 3 
work J 

vertical distance separating displacer from 
work piston at rest m 
peak allowable amplitude (semi-stroke) of 
displacer m 

peak allowable amplitude (semi-stroke) of 
work piston m 

vertical excursion of displacer from piston 

datum (rest) position m 

dimensionless volume—usually dead 

volume— V d /V sw - 

coefficient of dynamic viscosity Pas 

dead space parameter 

( = 8 r \nN T /(N T - 1) - 

density kg/m 3 

dimensionless mass, m/M , dimensionless 
mass rate, m'/coM 

dimensionless temperature, T/T c , e.g. T g = 
T g /T c - 

angular frequency rad/s 
regenerator volumetric porosity (pore 
volume/envelope volume) - 


mass of working fluid being cycled repeatedly through 
the same thermodynamic sequence; regenerative on the 
strength of the key, central component, the thermal 
regenerator. Comprising nothing more than a porous, 
conducting matrix with a high ratio of wetted area to 
volume, this passive, deceptively simple component 
stores and releases heat internally which would other¬ 
wise be exchanged with the surroundings at temperatures 
intermediate to the cycle limits of T E and T c . To this 
extent, the closed regenerative cycle fulfils the basic con¬ 
dition for thermodynamic reversibility. 

The cycle of operations is illustrated in Fig. 1, which 
shows in schematic form a cooler of coaxial configur¬ 
ation. A work piston controls total volume, while dis¬ 
placer and work piston in conjunction determine the dis¬ 
tribution of the working fluid between expansion and 



Subscripts 


ag 

air gap 

buff 

buffer space or ‘crank’-case 

E,C 

Expansion (temperature), Compression 


(temperature) 

c,e,r 

compression, expansion, regenerator 

d 

dead (unswept) volume displacer 

f 

(due to) friction 

ind 

indicated (as opposed to brake) 

P 

piston 

q 

(due to) heat transfer 

ref 

(constant) reference value 

rms 

root mean square 

rod 

relating to displacer drive rod 

sol 

solenoid 

sw 

swept (volume) 

w 

wetted (area) 


1. Introduction 




The Stirling cycle cooler operates on a closed, regen¬ 
erative cycle: closed because there are no valves, a fixed 


Fig. 1. Schematic representation of the sequence of volume changes 
of the Stirling cycle. The sequence is common between engine (prime 
mover) and heat pump (refrigerator) modes. 
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compression spaces. The regenerator, frequently com¬ 
prising a stack of fine wire gauzes, fills the hollow dis¬ 
placer and forms the flow path between the two variable- 
volume spaces. 

A typical cycle may be thought of as starting with the 
piston at the outer dead centre, and the displacer at the 
uppermost limit of its travel. The work piston rises, 
bringing about a compression of the entire working fluid 
mass, a rise in temperature overall and a corresponding 
overall rejection of heat. The majority of the working 
fluid mass being in the compression space during this 
process, the majority of the heat rejection occurs at or 
near compression space temperature, T c . The displacer 
is now caused to fall, and in the process brings about a 
transfer of a large fraction of the working fluid via the 
regenerator to the expansion space. The passing fluid 
deposits much of its heat in the regenerator, emerging 
close to expansion temperature, T B . Both piston and dis¬ 
placer now move downwards. The overall expansion 
lowers the temperature of the fluid occupying the expan¬ 
sion space. For as long as this temperature remains 
below that of the expansion space walls, the fluid 
absorbs heat, bringing about the required cooling effect. 
Towards the end of the expansion phase, the displacer 
is caused to rise, transferring working fluid via the regen¬ 
erator back to the compression space. In the case of an 
appropriately designed regenerator, the fluid picks up on 
this reverse pass the majority of the heat deposited in the 
matrix on the previous (hot-cold) pass, emerging with a 
temperature close to T c . 

The reversed Stirling-cycle cooler is capable of 
operating over a range of temperature differentials, and 
has been produced commercially in a wide range of heat- 
lift capacities. Miniature variants (Fig. 2) tend to operate 
between 77 K and ambient, lifting a fraction of a watt 
for an electrical input of a few watts. This impressive 
coefficient of performance has been achieved largely by 
driving the work piston by linear electro-magnetic exci¬ 
tation, thereby eliminating side thrust and resulting fric¬ 
tion. The displacer is suspended on a mechanical spring 
of some type and actuated by the pressure differences 
arising during the working cycle. 

The cycle of work, heat exchange and fluid flow vari¬ 
ations is a complex dynamic and thermodynamic interac¬ 
tion. Design of coolers from first principles calls for this 
interaction to be understood—a challenge traditionally 
approached by linearising the various disturbing and 
restoring forces on piston and displacer and treating the 
whole as a 2-degrees-of-freedom system. In the hands 
of de Monte and Benvenuto, the approach has been taken 
to a high art form [1]. 

On the other hand, it can be argued that most of the 
relevant phenomena are, in reality, non-linear. While lin¬ 
earisation of an isolated, non-linear phenomenon may 
frequently be justified, no-one knows for certain how 
good a representation of reality is achieved by arbitrary 


Detector substrate Variable-volume 



Fig. 2. Slightly simplified cross-section through 'B-Type’ cooler. 
Loudspeaker suspension is indicated for both piston and displacer. To 
gain full benefit from this approach to suspension and guidance it 
would probably be necessary to specify a ‘bucket’ piston suspended 
at both ends. Reproduced with the permission of DERA, Malvern. 

linearisation of the interaction of many non-linear 
effects. But perhaps the strongest case for a return to the 
basics of Stirling cooler dynamics is that the mathemat¬ 
ical contortions to which linear algebra leads assume 
greater complexity than a simple head-on solution of the 
original, non-linear problem, and the process involves an 
uncomfortably high ratio of algebra to insight. 

This paper accordingly adapts a well-known cycle 
analysis [2] to the case. Equations of motion for piston 
and displacer are written in terms of instantaneous gas, 
spring and drive solenoid forces. An implicit numerical 
integration scheme follows the gas processes and piston 
motions from start-up to cyclic equilibrium, revealing 
aspects of dynamic behaviour suppressed by the linear 
approach. Basic to the re-formulation of the problem are 
the equations of motion of displacer and work piston. 

2. Mechanical equations of motion 

The relevant features of the mechanical system are 
represented with associated notation in Fig. 3(a). Piston 
motion, y p , refers to the upper piston face, and the datum 
is piston rest position. The same datum is used for dis¬ 
placer motion, y d , which refers to the lowermost face 
of the regenerator housing. Amplitudes (semi-strokes) of 
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Fig. 3. Basic notation, (a) Schematic representation of gas circuit 
elements and related notation, (b) Coordinates for regenerator flow 
model. 


both members depend on operating conditions, but a 
notional maximum, Y amPp , is imposed on the piston by 
the step change in diameter. 

The equation of motion of the displacer takes account 
of two components of pressure difference: that between 
compression and buffer spaces (which acts on the cross- 
sectional area of the displacer guide rod, A xrod , and is the 
primary disturbing force), and that due to flow resistance 
through the regenerator. The latter acts on the whole dis¬ 
placer cross-sectional area, A xe . 

^dVd Pc(A X e A xrod ) /?(A xe T PbuffA xm d (1) 

- k d (y d - T sep ) 

The piston is subject to the difference in pressure 
between compression space and buffer space, to the 
restoring force of spring k p , and to the periodic driving 
force, F sol , of the solenoid. The latter is evaluated in 
Appendix A as a function of driving voltage and electro¬ 
mechanical characteristics. 

m p y p " = ip buff - /T)(A xc - A xrod ) - k p y p + F sol (2) 

In Eqs. (1) and (2) the double prime (") signifies the 
second differential with respect to time, t. Instantaneous 
compression space volume, v c (<f>) is given by 


v c ((p) = (T amPp - y p )(A xc - A xe ) (3) 

A (y d y p )(A xe A xrod ) 

With K ampd for maximum y d before the displacer contacts 
its stops (or the spring becomes coil-bound) instan¬ 
taneous expansion space volume is 

= (T amPd - >’ d )A xe (4) 

Instantaneous buffer-space pressure, p buff , is estimated 
by assuming adiabatic compression and expansion of the 
volume of fluid below the piston. With V 5uff for the value 
of the volume with piston and displacer at rest and 
recalling that positive y p and y d bring about an increase 
in volume: 

_ _ ^hulT _ _ 

Pbu,f P " f I V bn „ + y r (A xc - A rod ) + 

The volume variations may be thought of as boundary 
conditions for the thermodynamic statement of the prob¬ 
lem—which in turn determines the pressures and hence 
the motions, y p and y d . 

The context calls for a tailored approach to cycle 
simulation: dynamics at start up (T F = T c = ambient) 
promise to be different from those at rated conditions— 
and, indeed, from those at intermediate temperatures. In 
real time, some thousands of cycles are involved in coo¬ 
ling from ambient to operating temperature, justifying 
some corner-cutting with temperature solutions provided 
reasonable fidelity is retained to flow conditions, since 
it is these which determine instantaneous pressure drops. 
The thermodynamic treatment proposed here builds on 
Finkelstein’s Generalised Thermodynamic Analysis [2], 
but avoids certain problems of the original formulation: 
preparing the defining equations as a matrix of 
unknowns rather than substituting and eliminating 
algebraically keeps the expressions for variable volumes 
out of the denominators, thus avoiding division by zero. 
Superimposing the steady-flow regenerator solution [3] 
at each integration step provides an approximation to 
regenerator thermal response which is sensitive to the 
parameters of the real component—hydraulic radius, r h , 
volumetric porosity, % etc. A technique is incorporated 
for calculating regenerator pressure drop simultaneously 
with the pressure and temperature solutions rather than 
imposing it retrospectively in terms of computed flow 
rates. 

3. Gas process model 

3.1. Energy equation for variable-volume spaces 

Either variable volume, V of Fig. 3(a) contains mass, 
m, instantaneously at temperature T g and pressure p. An 
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increment of time d t is considered corresponding to RT g dm + Rm dT g — V dp = p dV (8) 

increment dc/> in dimensionless time (or ‘crank’ angle). 

Energy conservation for the space may be expressed as The dV term is on the right-hand side, reflecting the fact 

that the dV are known via Eqs. (3) and (4) from the 
dQ + dH - dW = dU equations of motion. 


In general, m changes by dm in the interval: if mass 
enters it does so from the regenerator and at the relevant 
regenerator outlet temperature, T ro . When dm leaves, on 
the other hand, it does so at the instantaneous value of 
the varying T g . Following Finkelstein [2], step function 
U{dm) is defined such that: 

U (d/?7,) = 1 for dm > 0 
= 0 for dm < 0 

This permits the effects of dm to be handled in a single 
algebraic term. With h for instantaneous coefficient of 
convective heat transfer and A w for the corresponding 
wetted area, the energy balance may be written as 

M W (T W - T g ) dt + c p dm{r g [l - U(dm)] (6) 

+ U(dm)T ro } — p dV = c v (T„ dm + m d T g ) 

3.2. Mass conservation 

On the basis of a linear variation of temperature 
between T E and T c and a small pressure gradient, the 
instantaneous mass content of the regenerator is m r = 
1/2 (p e + p c )V dI In N t /{RT c (N t — 1)}[2,3]. Mass inven¬ 
tory, M, of the gas circuit being assumed constant, the 
sum m e + m c + m r = M, may be notionally differen¬ 
tiated to give 

V dr In N t 

dm e + dm c + ^ l/2(d p e + dp c ) = 0 (7) 


3.3. Gas law 


3.4. Pressure drop in terms of flow resistance 

The steady-flow component contributes the majority 
of the pressure drop across the regenerator. A severe 
temperature gradient and a substantial phase shift 
between the flows at the two extremities suggest inte¬ 
grating local pressure drops over elemental lengths, Ax, 
of regenerator packing to form instantaneous overall 
pressure difference. To achieve this in terms of the vari¬ 
ables of the implicit solution, namely, the two pressures, 
p e and p c , and the corresponding mass rates, mj and m c ', 
calls for a bit of sleight of hand. 

The usual expression for steady-flow pressure drop 
takes a negative sign (flow down the pressure gradient). 
The present treatment requires p e — p c to be positive 
when flow is from expansion space to compression. With 
reference to Fig. 3(b): 


doc u 

Pj+i ~ Pj = + l/2pw C f — n 

r h \u\ 

The term u/\u\ keeps the higher pressure upstream of the 
lower. The expression may be re-written in terms of 
local m' = puA ff . With p = p/RT and with x/L reg = A 


m' 2 C f L ieg RT u 

Pe + Pc Viff | U 


dA 


At any given time, m' varies with location, x, between 
mf and m c ', and over substantial intervals of crank angle 
may be inwards from both ends of the regenerator or 
outwards from both ends. An expression for instan¬ 
taneous m' at location A which accommodates all such 
flow possibilities but with some compromise at locations 
intermediate to the two ends is: 


An advantage of formulating in terms of a matrix of 
unknowns is that it is unnecessary to make algebraic 
substitutions in terms of working fluid properties from 
the gas law. Fluid temperature falls below nominal T E 
during expansion so that, depending on the gas, con¬ 
ditions may become sufficiently close to the critical point 
to warrant consideration of an equation such as that of 
van der Waals. 

For convenience of presentation the treatment pro¬ 
ceeds in terms of the ideal gas equation, but will be seen 
to be capable of ready adaptation to more sophisti¬ 
cated models. 

The gas law (pV = mRT) is differentiated, the result 
applying to either working space: 


m'{ A) = — m' e + {in' c + m' e )A (9) 

m' e ( — 1 + A) + m' c A 

Eq. (9) reflects the fact that positive mj is flow into the 
expansion space, and therefore flow out of the expansion 
end of the regenerator. 

The customary approximation to local temperature 
is linear: 

T{ A) = T e + (T c + T e ) A (10) 

Taking out constant factors, noting that m' is dm/dt, and 
evaluating pressure drop at mid-interval dr: 
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Pe + 1/2 d Pj + 1 - p c - 1/2 dpj = 

L reg R 1 (dm e ( — 1 + A) + dm c A) dm 

r h Pe + p c Af f 2 dt 2 | dm| 

...C f (r E + (T c + r E )A} dA 

Summing over the length of the regenerator gives an 
expression which may be made quasi-explicit in the dm e 
and dm c as well as in the pressure increments: 

L reg R 1 dA 

1/2 d pc ~ 1/2 d pc = ——- — — {dmddmj:, 

r h Pe + Pc A ff - d/ 1 " 

+ dm c iy... + dm c [dm e Z 2 + dm c Z 3 ]} - p c + p c 

n 

E. = 2 ( - 1 + A) 2 C,(r E + (T c + r E )A} dm/|dm| 

j = 1 

n 

E 2 = 2 A 2 C f {7- E + (T c + r E )A) dm/|dm| (11) 

j = 1 

n 

E 3 = 'Z A( - 1 + A)C f {T E + (7 C + r E )A) dm/|dm| 

.7 = 1 

In evaluating the L the C f and dm/|dm| corresponding to 
local, instantaneous m’ and T (Eqs. (9) and (10)) are used 
to form the local, instantaneous Reynolds number, N re . 
The corresponding friction factor, C f , then follows from 
an interpetation of the Kays and London correlation [7] 
for the appropriate value of volumetric porosity, %. For 
the % of the specification: 

C f = 40/iV re + 0.25 (12) 

3.5. Regenerator temperature distribution 

Until 5 years ago, temperature solutions for regener¬ 
ator response to the flow and pressure conditions of the 
Stirling machine were unavailable. This has changed 
since development of a ‘natural coordinate’ system [3] 
for creation of the coincident Lagrange/Euler integration 
grids required for simultaneous solution of fluid and 
matrix temperatures. Acquiring a comprehensive picture 
of regenerator in situ transient response for the crank- 
driven Stirling cycle machine is now a routine matter. 
On the other hand, it calls for generation of a fluid par¬ 
ticle trajectory map for each set of operating conditions. 
While the computational effort was not great for the 
original formulation [3], the generalisation of ‘natural’ 
coordinates called for by the arbitrary volume and tem¬ 
perature variations of the present treatment [4] require 
significantly more numerical processing—hence the 
search for a suitable simplified approach. 

Thermal capacity ratio, /V TCR is defined in terms of 
wire properties, p w , c w as p w c w T vet /p ref . The numerical 


value in the present instance is 1250. This is twice that 
of the Philips MP1002CA air engine, which in turn 
exceeds that of the GPU-3 engine by a factor of almost 
5. This, in conjunction with the massive volume of 
regenerator material in the cooler relative to that of the 
prime movers, suggests provisionally ignoring matrix 
temperature swing as insignificant in relation to fluid 
temperature fluctuations. Actual values are computed 
later. 

With matrix temperature invariant, the variation of 
fluid/wall temperature difference, At = (77 — T w )/T ref , 
along a particle path is given [3] by 

^ + P(<f>) At = (13) 

In Eq. (13) 

p = |w|/v st — = |m|ntu 

r h 

The NTU are evaluated via N st corresponding to the local 
instantaneous Reynolds number, N re , using the Kays and 
London [7] correlation appropriate to the value of volu¬ 
metric porosity, % 

N„N p , m = a/Nj’ (14) 

Q takes account of the pressure-swing term: 

_ 5t w t w y ~ 1 gjr 

5A ip y dcp 

The analysis is presented in full in [3]. For an interval 
over which flow may be considered steady, the solution 
of Eq. (13) proceeds in steps starting at the inflow end. 
The first step sets At to inlet temperature difference: 

At^ + A4) = Q/P + (At^ - Q/Pff~ p ^ (15) 

Inlet temperature difference, AT in , is established, by 
definition, at a time previous to that of the current com¬ 
putational step, and so is always known regardless of 
flow direction. The only subtlety in the computation is 
that it is necessary to ‘step back’ by the appropriate 
angular amount of dimensionless time, cjj, for the AT in 
appropriate to calculation of the At at given location A. 

The solution for A = 1 is compression space outlet 
temperature difference, that for A = 0 the expansion-end 
outlet difference. These provide the working space inlet 
temperatures, T ro , required for evaluation of Eq. (6). Eq. 
(15) also provides the At intermediate to the two ends 
for plotting the temperature envelopes and reliefs dis¬ 
cussed later. The NTU on which the P are based are 
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calculated in terms of the Reynolds number, N re , corre¬ 
sponding to local instantaneous flow conditions. 

The approximation improves as regenerator dead 
space ratio reduces and flow at outlet falls increasingly 
in phase with that at inlet, i.e. as steady flow conditions 
are approached. 

3.6. Instantaneous heat transfer coefficient in the 
variable-volume spaces 

Values are required of the coefficient of convective 
heat transfer, h , for use in the energy equations (Eq. (6)). 
The task is best done in terms of Stanton number, N st , 
correlated in turn with Reynolds number, N re (and 
Prandtl number, V pr ) for conditions in the variable-vol¬ 
ume spaces. The only resource known to the writer is 
Annand’s correlations [5] for conditions in the cylinder 
of the reciprocating internal combustion engine. Annand 
determined values of a and b for the following 
expression for heat rate, q 

q' = A W '\a | R b e (T w - T g ) + c(K - T$\ 


working spaces is coS normalised by iro>L ref , or ShrL ref , 
where S is stroke—of displacer or piston as appropriate. 


4. Preparation for solution 

When Eqs. (6) and (8) are written for each of the two 
working spaces, the resulting total of six equations are, 
but for the U terms, linear in the six unknowns, dm e , d m c , 
d p e , dT ge , dT gc , and d p c (in that order). For generality, 
computation proceeds in terms of normalised variables 
and dimensionless parameters. Appendix B carries out 
the normalisation and lists the coefficients of the 
unknowns. 

Individual solutions (increments in the dependent 
variables) are added to respective values of <x e , <x c , i// e , 
Tge, T gc and i f/ c (in that sequence) at each integration step 
to form values corresponding to new crank angle f + 
df. The formulation is free of division by terms which 
can take a value of zero, and consequently there are no 
problems with singularities. When f = 2ir, values are 
compared with those for f = 0, as many cycles being 
repeated as are necessary for convergence. 


in which the Reynolds number, R e , is defined as 
R c = IfSpDIp 


5. Discretisation and normalisation of motion 
equations 


Neglecting the radiation component as being insignifi¬ 
cant at the temperatures involved, and defining h as the 
ratio of heat rate, q', to the product of wetted area, A w , 
with temperature difference, T w — T g , leads to 


h = 


akRf 

D 


(16) 


Formulation of the numerical solution proceeds by 
forming a finite-difference approximation to the differ¬ 
ential coefficients of the motion equations. For example 


d 2 y p _ Ay p ew - Ay° ld 
d f~ A t 2 


(18) 


Normalisation of the energy equation (Eq. (6)) yields a 
term hAJcoMc p , which is recognizably NTU—although 
not yet in useful form. However, multiplying numerator 
and denominator by \u\L ref mco brings about a conversion 
to NTUcr|u|, where there is now the ‘proper’ NTU = 
N st L/r h . It remains to define hydraulic radius, r h , for the 
working spaces. 

r h is defined as wetted volume, V w , divided by wetted 
area, A w . When exposed cylinder wall height is h, V w 
for the expansion space is l/47iD 2 /r, while corresponding 
A w is 2{ 1/47 tD 2 } + TiDh. Forming the ratio, simplifying 
and normalising: 

/V 1 = 2/fi/L ref + 4/D/L ref (17) 

r h for the compression space is dealt with in a similar 
way. For consistency with Annand’s definition of Reyn¬ 
olds number, dimensionless velocity, u, required for both 


In this expression the increment of motion Ay p new is the 
‘unknown’ to be determined by solution of the complete 
set of simultaneous difference relationships. Piston 
dynamic response may be thought of as being centred 
on current time, t (or dimensionless time or crank angle, 
f = cot) as suggested by the inset in Fig. 3(a): the forces 
responsible for motion (spring, pressure difference, 
solenoid drive force) are those representative of the 
interval, i.e. mid-interval values. Given that solution is 
proceeding in the sequence / — 1 to unknown i + 1 via 
i, mid-interval values are ‘known’. This leaves Ay p new as 
the sole unknown in the motion equation, which is there¬ 
fore explicit in that unknown, and does not require to 
be part of the matrix of the unknown thermodynamic 
increments, dm e , d m c etc. 

The same normalising parameters are applied to the 
motion equations as to the gas process model, and the 
unknown in each case made explicit. With A p for dimen¬ 
sionless piston excursion, y p /L ref 
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dA p 


new 


dA p old + 


A cf) 2 
(*V<) 2 



+ 


PrefL, 


'ref 


#p 


Ob 


*Ac)(a xc 



Table 1 

(19) Selection from the data forming specification and operating point 
which produced Figs. 4, 5, 7-9. All items are believed to be achievable 
in the cooler illustrated in Fig. 2. Flowever, with the exception of exci¬ 
tation frequency, /, which is cited with the permission of DERA Mal¬ 
vern, all ‘data’ are assumed, and do not represent the official specifi¬ 
cation of any cooler 


In Eq. (19) F' ol is an abbreviation for the mid-interval 
(normalised) solenoid driving force derived in Appendix 
A. The corresponding expression for the increment in 
displacer motion is 


dA3 ew = dA° ld + 


A cjr 


— A d + 


PrefL 


ref 


(o)/o ) np ) 2 [ “ a # p 


( 20 ) 


[(*Ac(Axc *-Arod) t/f e rt X( j "f t/ffAxTod)] 


Evaluating Eqs. (19) and (20) in advance of solving the 
6X6 matrix of ‘thermodynamic’ unknowns provides 
values for the right-hand sides (volume changes) of Eq. 
(4)a and b by simple differentiation of Eqs. (3) and (4) 
for the respective volumes. 


6. Specimen simulated performance 

The writer has permission for the reproduction of Fig. 
2 of the DERA B-type cooler, whose excitation fre¬ 
quency of 72 Hz is unclassified information. Pending 
specific permission to disclose dimensions, masses, etc., 
the entries in Table 1 are consistent with, and believed 
to be achievable in, the cooler as shown, but not, except 
by coincidence, part of the official specification. A 
woven-wire mesh regenerator is assumed. A wide range 
of relevant gauzes is manufactured by Messrs. G. Bopp 
and Co. [6] and can be supplied punched to size. 
Assuming rectangular, square weave justifies use of fric¬ 
tion factor-Reynolds number (C t ~N re ) correlations from 
the Kays and London compendium [7]. There is no 
attempt to ‘grade’ the packing despite the marked vari¬ 
ation in physical properties and flow conditions between 
hot and cold ends, although the option has afforded the 
basis of an interesting study [8]. 

The cooler itself is the result of a programme of 
experimental development work unconnected with the 
current theoretical study. It is therefore of interest to see 
what, if anything, the latter has to say about the former. 

The cooler is ‘started up’ with expansion and com¬ 
pression space temperatures, T B and T c pre-set to ‘rated’ 
values (Table 1). The starting condition for the regener¬ 
ator is a distribution of matrix and fluid temperature 
which is linear between T E and T c . Fig. 4 shows indi¬ 
cator diagrams: the smallest trace is that for the expan¬ 
sion space, the largest for the compression space. The 
rightmost loop is the net indicator diagram. 


Working fluid 

helium 


p u . f charge press 

8 

bar 

B 

2.0 

T(esla) 

f excitation freq. 

72* 

cps 

k d 

5218 

N/m 

k p 

841 

N/m 

m d 

8.0E - 03 

kg 

m p 

25E - 03 

kg 

N 

125 

solenoid turns 

T c 

300 

K 

t e 

77 

K 

v sw 

2.13 

cm 3 

Derived data 

A re f( = O 

1.65E - 04 

m 2 

8e( = F de /V sw ) 

0.05 

(nom) 

s r ( = VJV S w ) 

0.447 


S c ( = E dc /V sw ) 

0.05 

(nom) 

L ret ( = O 

0.0128 

m 

L ie f/L ref 

4.8 

- 

^MA 

0.0073 

- 

Nsg 

1.04E + 08 

- 

VAef 

0.003 

- 

n t = t e /t b 

0.256 

- 

'no-: 

1250 

- 

ftj/tonp 

1.0 

- 

co/M nd 

1.4 

- 

% 

0.686 

- 



O.o 1.0 v/F sw 2.0 

Fig. 4. Specimen computed indicator diagrams. The smallest, left¬ 
most is for the expansion space. The larger trace at the left is for the 
compression space. The rightmost loop is the overall indicator dia¬ 
gram, and represents indicated work per cycle. 


Fig. 5 shows computed variation of working space 
temperatures, T ge and T gc over a cycle (0-2 tt rad). Rela¬ 
tively small swings either side of nominal cycle tempera¬ 
ture limits confirm that high values of NTU e and NTU C 
are being achieved at ‘rated’ operating frequency and 
charge pressure chosen. The swings in temperature are 
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Fig. 5. Cyclic variation of dimensionless working space tempera¬ 
tures, T ge /T c and T g JT c . Note the very small swings, consistent with 
the high rms values of ‘number of transfer units’ (NTU) for the 
exchange processes in the variable-volume spaces. 

in phase with those of pressure (dashed curve), reflecting 
an element of adiabatic heating. The lower (expansion) 
trace is biased below nominal expansion-space tempera¬ 
ture of N t , consistent with heat flow to the working fluid. 
The upper trace is biased above, consistent with heat 
rejection at sink temperature. 

Computed regenerator NTU are so high that the fluid 
temperature envelope plots as a straight line between 
cycle limits. To confirm that the temperature solutions 
reflect the picture given by more comprehensive treat¬ 
ments (e.g. those of [3]) NTU have been arbitrarily 
reduced by a factor of 50 for Fig. 6, in which (a) is a 
set of superimposed fluid temperature envelopes for a 
complete cycle at cyclic equilibrium, and (b) is the corre¬ 
sponding fluid temperature relief. 

Fig. 7 shows normalised mass flow rate in function of 
‘crank’ angle at opposite ends of the regenerator. The 
larger amplitude of the expansion trace is consistent with 
the greater fluid density, but neither trace bears much 
likeness to mass-rate records which have been computed 
for Stirling prime-movers, such as those of Creswick [9]. 

Fig. 8 plots the cyclic variation of Reynolds number, 
AU, at the regenerator extremities. The N re are corrected 
for the temperature dependence of coefficient of 
dynamic viscosity, /r, by use of the Sutherland Law [3]. 
The substantially larger values in the expansion space 
reflect high density (in the numerator of N Te ) and reduced 
coefficient of dynamic viscosity (in the denominator). 

In Fig. 9 the lowermost, light trace is driving voltage, 
which leads piston motion (lower heavy trace) by about 
45°. The central light trace is pressure difference. The 
upper light trace plots expansion space pressure vari¬ 
ation. The upper heavy trace is the lowermost face of 
the displacer. At the operating condition modelled, the 
amplitude of the regenerator pressure drop is about one 
third of expansion space pressure excursion. In a Stirling 
prime mover this ratio would be excessive and would 
probably prevent operation anyway. In the cooler as 
modelled it represents the cyclic interaction necessary 


(a) 



Fig. 6. Fluid temperatures computed from simplified regenerator 
theory for comparison with comprehensive solutions (see, e.g. [3]). (a) 
Temperature envelopes for working fluid occupying the regenerator, 
exaggerated by decreasing local, instantaneous NTU by a factor of 50. 
The machine fails to lift heat with the regenerator ‘modified’ in this 
way—as reflected by the fact that the temperature envelope lies perma¬ 
nently above the lower limit of the cycle, i.e. above N T (heat exchange 
within this component is so effective that to obtain this plot it was 
necessary to decrease computed NTU for each integration step by fac¬ 
tor of 50). (b) Fluid temperature relief corresponding to (a) above, i.e. 
NTU decreased by a factor of 50 to exaggerate temperature excursions. 



Fig. 7. Computed variations of dimensionless mass rate, <x e and cr c . 
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Fig. 8. Computed variations of Reynolds numbers at regenerator 
extremities. 



Fig. 9. Specimen motion diagram. Lowermost curve is driving volt¬ 
age shown leading piston motion by about 45°. The bold trace above 
this represents the motion of the uppermost face of the piston. The 
light superimposed trace is pressure difference, Ap = p e - p c . The 
upper heavy trace follows the lowermost face of the displacer. The 
uppermost light trace shows expansion space pressure swing. Note that 
pressure difference is a sizeable fraction of pressure swing. 


for achievement of the required phase between volume 
variations. At the same time it suggests where a search 
might be directed for a less inefficient way of realising 
the Stirling cycle in the free-piston embodiment. 

7. Deductions 

In overall terms, the computer models behaves in 
accordance with the supposed mechanism of operation: 
for example, increase of displacer drive rod diameter 
causes the displacer to over-travel, impacting the stop at 
the upper limit and/or the piston in the lower, as in Fig. 


10(a). Excessive drive voltage, E 0 , on the other hand, 
causes over-travel of the work piston (Fig. 10(b)). In 
both cases indicator diagrams become distorted, and 
traces of mass rate and Reynolds number reflect the sud¬ 
den changes in flow direction. Raising flux density, B, 
across the airgap brings about a rise in electro-mechan¬ 
ical conversion efficiency, and calls for a higher driv¬ 
ing voltage. 

The model has not yet been subject to mechanised 
optimisation, but after a large number of manual pertur¬ 
bations of the variables, operation at the parameter 
values of Table 1 appears to be close to optimum (or to 
an optimum): computed heat lift at 77 K is 1.73 W for 


(a) 






N. / 













\ / \ / \ / 



Fig. 10. Features of ‘off-design’ performance. Any mechanical con¬ 
tact appears to reduce heat lift capability and electro-mechanical con¬ 
version efficiency, (a) Displacer ‘over-driven’ by excessive value of 
d mA . (b) Piston ‘over-driven’ by excessive input voltage. 
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a mechanical input of 29.37 W—or 16.97 WAV lifted, 
corresponding to a thermo-mechanical COP of 5.89%. 
Hymatic Enginering Co. claim [10] to lift 0.373 W from 
82.2 K for 10.77 W (electrical) input, or 28.87 WAV nett 
in a somewhat smaller cooler (nominal 7 mm diameter 
cold-finger). To compete, the cooler of the present study 
would have to achieve an electro-mechanical conversion 
efficiency of 10.77/28.87 = 37.3%. In fact, the 
efficiency computed on the assumption of the 2.0 T(esla) 
flux density through the solenoid is only 22%, remedi¬ 
able only by an improved permanent magnet specifi¬ 
cation, and/or ironwork giving higher flux concentration. 

The study affords less obvious insights into design 
and operation: 

• A frequently-expressed aim of work to develop this 
type of device is enhancement of regenerator thermal 
performance. The finding here is that temperature 
recovery ratio is not an issue: the temperature differ¬ 
ence between fluid stream and matrix attains its equi¬ 
librium value of 0.18°C rms within the first inte¬ 
gration step of the computational sequence through 
the regenerator. Step size used here is standardised at 
1/20 of regenerator length, so regenerator heat transfer 
capability is more than adequate, as confirmed by a 
high rms value of NTU r = 132 at the rated conditions 
assumed. The massive value of thermal capacity ratio, 
/V TC r (1250 from Table 1) contributes to a high tem¬ 
perature recovery ratio. 

• The expectation from experience with Stirling prime 
movers is for heat transfer in the working spaces to 
be poor. By contrast, rms temperature difference in 
the expansion space is 6.0°C; that in the compression 
space 13.5°C. Had this information been available at 
the time the NAX-106 cooler was designed [3], the 
compression-end heat exchanger would not have been 
incorporated, and dead space, compression ratio and 
overall performance would have benefited. 

• The basic thermodynamic design problem is achieve¬ 
ment of a displacer motion lagging that of the piston 
by a suitable amount. All other things being equal, 
displacer amplitude is controlled by drive rod diam¬ 
eter—with changes in buffer space dead volume also 
having a small effect—but the phase relationship is 
a result of a more subtle combination of factors. To 
appreciate these it is necessary to be clear that for 
present purposes the natural frequencies, cn p and u> d , 
of piston and displacer are those calculated in terms of 
respective masses, m p and m d , and the corresponding 
stiffnesses, k p and k d of the mechanical springs and 
not the pseudo-spring constants of the gas circuit of 
the ‘linearised’ approach. 

Thus it is seen that increasing the mass of the dis¬ 
placer (say) whilst maintaining natural frequency (by 
increasing k d in proportion) retards response to the 
pressure difference acting on the drive rod and 


increases phase shift. 

A decrease in hydraulic radius, r h , of the wire mesh 
may be used to increase phase shift, but generally at 
a cost of increased pumping power. In a conventional 
Stirling engine or cooler, i.e. one with kinematically- 
synchronised piston and displacer, the increase would 
be thought of as an unacceptable performance penalty. 
A different interpretation applies here: the traces show 
a pressure drop to be a sizeable fraction of pressure 
swing. If this is the price to be paid for achieving a 
phase shift which in turn achieves heat lift, then it has 
to be paid. On the other hand, with the generous over¬ 
supply of regenerator heat transfer capacity suggesting 
an increase in hydraulic radius and corresponding 
reduction in pumping power, a search for performance 
improvement might focus on other parameters having 
an influence on phase shift. 

• Electro-mechanical conversion efficiency is a strong 
function of operating conditions, and for the specifi¬ 
cation of Table 1 peaks at about 25%. The most 
influence parameter of conversion efficiency is flux 
density, B. The value assumed here is 2.0 (from the 
demagnetisation curve of neodymium-iron-boron) 
enhanced by a factor of 1.53 for the concentrating 
effect of the ironwork. Use of a higher value calls for 
increased driving voltage, E 0 , but considerably 
enhances conversion efficiency. Where overall coef¬ 
ficient of performance (ratio of heat lifted to power 
input) is critical it appears important to embody a 
magnet material of high remanence and/or a high ratio 
of magnet cross-section to cross-section of the iron¬ 
work. 

• For many combinations of design parameters it is 
possible to determine by trial and error a combination 
of input voltage, E 0 , and angular frequency, co ( = 
2nf) which will force a favourable phase-shift 

between piston and displacer, and a corresponding lift 
of heat from the expansion end. If development is 
approached in this way, the resulting heat lift rate (W) 
and overall COP are not under the control of the 
designer to the extent that they might be. 

Availability of a computer simulation, by contrast, 
allows selection of the desired combination of design 
parameters and operating conditions in advance of 
manufacture. 


Appendix A 

Elementary view of the electro-magnetic circuit 
Supplementary notation 

B magnetic flux density = p qH in Tesla 

free space 

D diameter of coil m 
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d diameter of wire of coil m 

E 0 amplitude (semi-swing) of supply V 

voltage 

e emf V 

e h back emf V 

i current A 

L ag length of air gap m 

L c overall physical length of coil m 

/ length of conductor m 

F instantaneous force N 

H magnetic field strength, or amp/m 


‘magnetising force’ 

N total number of turns in solenoid 

n number of turns cut by flux path 

/x 0 permeability of free space 477 X Henry/m 

1CT 7 Henry/m 

[jl v relative permeability—ratio of 

permeability of given material to 
that of free space: B = p^^qH 
p resistivity ohm.m 

co angular frequency rad/s 


flux path energised by a permanent magnet and a coil 
wound into a rigid annulus. Permanent magnet and flux 
path are in the form of continuous annuli. The coil is 
rigidly attached to the work piston which is free to oscil¬ 
late along the axis of the assembly. The axial length of 
the coil, L c , exceeds that of the pole piece by an amount 
which ensures that the flux path cuts the same number 
of turns regardless of coil axial location. This condition 
confines coil excursions to the limits indicated. 

Motion of the coil at velocity u, generates a back 
emf, e b : 

e b = t tBuDNW/L c (Al) 

At the sort of cyclic frequencies envisaged, the corre¬ 
sponding instantaneous current is given by the ratio of 
net voltage to total coil resistance: 



(A2) 


Subscripts 
ag air gap 

b back (emf) 

c coil 

sol solenoid 

w wire 


Principle of operation 


Instantaneous axial force generated is proportional to 
current (F = Bil) 

F sol = jrBiDNW/L c (A3) 

It is necessary to have E sol normalised to the same base 
as the equations of motion, i.e. by k p L ref = k p V sw 1/3 . 
Assuming a supply voltage of the form e = E Q sin cot, 
substituting Eq. (A2) into Eq. (A3) and normalising: 


Fig. 11(a) is a schematic representation of an electro¬ 
magnetic circuit. This comprises a stationary magnetic 


F so i 

£pKw 1/3 


ttBDNW 1 

kpVlwLc R 


{E 0 sin cot 


irBuDn W/F c } 





Fig. 11. Notation for electromagnetic circuit, (a) Solenoid with n 
effective turns cut by flux density B. (b) Part of demagnetisation curve 
of permanent magnet. 


There is some economy of algebra in expressing coil 
resistance, R, in terms of resistivity of wire material and 
coil geometry: 


irnDp 
1/4 m/ 2 


ANDp/d 1 


(A4) 


Making the substitution, and normalising u as per the 
thermodynamic treatment, i.e. by u = u/coL re{ = 
dA p /d 4>, where dA p is the current, instantaneous linear 
increment in piston position 


^sol 

k p Vl£ 


DGl{sin (jo — DG2 u} — F' sol 


E 0 BWird 2 

Apk p VlL 3 L c 


DG2 


TrBooL ret DNW 

E 0 F C 


(A5) 
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Determination of operating point 

Fig. 11(b) represents part of the demagnetisation curve 
(. B-H) curve of a permanent magnet. Over the range for 
which the curve is linear, hysteresis may be ignored, i.e. 
a shift in operating point, Q, caused by the demagnet¬ 
ising effect (air gap, ampere turns or an opposing mag¬ 
netic held) may be reversed by removal of the demagnet¬ 
ising effect. Point Q is determined in terms of the B and 
H allowed by the hux circuit of which the permanent 
magnet is part. 

Law of flux continuity 

In Fig. 12(a) the flux, 4>, which links the permanent 
magnet passes through the soft-iron circuit and the air 

gap 

^ _ ^ ma gnet^x m agnet _ ^iroi/fx imn — ^aiiAx air (A6) 

Ampere’s circuital law, jH-dl = ni 
Applying the law to Fig. 12(b) 


Eq. (A8) plots into Fig. 11(b) as a straight line of nega¬ 
tive gradient, determining operating point Q. 


Appendix B 

Preparation of gas process equations for solution 

When Eqs. (6) and (8) are written for each of the two 
working spaces, the resulting total of six equations are, 
but for the U terms, linear in the six unknowns, dm e , 
d m c , d p e , dF ge , d T gc and d p c (in that order). 

For generality, computation proceeds in terms of nor¬ 
malised variables and dimensionless parameters. Eq. (2) 
is written out for each working space in turn: 

NTU e (<T J8 e )(N T - T ge ) d<p + do- e {T ge [l 
- f/(dcr e )] + U(da £ )N T } - ip e d8 e - T ge dcr e (Bla) 
+ cr e dT ge = 0 

NTU c (o- c /fl c )(l - T gc )d(f) + dcr c { T gc [ 1 

U(da c )] + f/(dcr c )} - ifj c dS c - T gc dcr c (Bib) 

+ cr c dT gc = 0 


iion^iroii "F iii]A a g HI (A7) 

In the annular iron path shown L jron = ttD — g ~ ttD. 
For the special case of ni = 0 (n = 0 and/or i = 0) 
consideration of remanent fields leads to H ilon 7rD = — 

ZJ T 

** airbag* 

In Fig. 12(b) there is no current-carrying coil, so ni 
0 and F/ magnet L magnet A H^ ron L^ on A H a ^ r L ag 0. Taking 
account of Eq. (A6), (A7), and using jU 0 )U r air = 
gives 


The NTU are respective ‘number of transfer units’, typi¬ 
fied by: 


NTU e 


aK ~ 

ttN pr L ref 


Normalised mass balance is 

dcr e + dcr c + r , D l/2(di// e + di// c ) = 1 (B2) 


B 


JH n 


magnet'' - 1 magnet 


f t o(A Xair /A Xmagnet ) (T magnet /L ag ) (A8) 


The normalised description is completed by writing the 
difference form of the gas law for each working space: 



(b) i/zLj ron 



Fig. 12. Determination of operating point in terms of circuit 
geometry: (a) ironwork circuit; (b) air gap and ironwork as magnetic 
circuit elements. 


T ge dtJe + cr e dT ge “ diff = if/ e dfi e (B3a) 

Tge d<r c + cr c dT gc ~ 8 C diff = if/ c d5 c (B3b) 

The integration process involves writing the coefficients 
of the six unknowns d<r e , dcr c , dif/ e , dT ge and dT g and dif/ e 
(in that sequence) as a 6 X 6 matrix and solving numeri¬ 
cally. Non-zero coefficients of unknowns in the order 
just given are: 

Energy-expansion end 

a( 1 , 1 ) = T ge ( 1 - U(da e ) - 1/y) + U(dcr e )T roe 

a(l,4) = - l/2NTU e Ac(> - ajy 

b( 1) = - NTU e (A T - T ge )Ac/) + ijj e dv e (y - l)/y 
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Energy-compression end 

a(2,2) = T gc (l - f/(d<x c ) - 1/y) + U(dcr c )T roc 

a( 2,5) = - l/2NTU,aA4> - ajy 

b(2) = - NTU C (1.0 - r gc )A <f> + ijj c dv c (y ~ l)/y 

The T roe and T roc are the appropriate regenerator outlet 
temperatures computed with the aid of Eq. (15) of the 
main text. 

Mass conservation 

a(3,l) = 1 

a( 3,2) = 1 


Momentum 

a(6,l) = - (da- e X x + da- c X 3 )DGM 
a(6,2) = - (da c X 2 + da- e X 3 )DGM 
a(6,3) = 1/2 
a(6,6) = - 1/2 
b( 6) = ~ ijj e + ijj c 

(\L re rj A 1 1 

DGM = r reg — — VV2 , , , 

L ref a ffr 2 r h A(/>- ifj e + <// c 

The various are the normalised equivalents of the %„ 
of Eq. (12). 


(3(3,3) = 1/21^, 
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